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Abstract 

We extend the mathematical model based on stochastic differential 
equations describing the error gained by an atomic clock to the cases 
of anomalous behavior including jumps and an increase of instability. 
We prove an exact iterative solution that can be useful for clock simu¬ 
lation, prediction, and interpretation, as well as for the understanding 
of the impact of clock error in the overall system in which clocks may 
be inserted as, for example, the Global Satellite Navigation Systems. 
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1 Introduction 

Atomic clocks are usually the heart of complex scientihc and technological 
systems and clock behavior has a direct impact on the overall performance of 
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the system. For the study, simulation, analysis, characterization, and inter¬ 
pretation of all these systems it is therefore mandatory to have a consistent 
mathematical clock model able to represent its typical behavior. Several 
studies have dealt with this matter P El El H mil [7]. 

Recently, particularly in space applications such as the Global Satellite 
Navigation Systems (GNSS), anomalous clock behaviors have been experi¬ 
mentally observed [HlIHlIiniElElEniE]: phase or frequency jumps, increase 
of instability, general non stationary trends. Such anomalous behaviors have 
an important impact on the global system and may affect the overall per¬ 
formance very seriously. This could have, in case of GNSS, dramatic conse¬ 
quences even security of life. 

In this paper we propose a mathematical model including the typical 
anomalies of atomic clocks and we show how the stochastic multidimen¬ 
sional differential equation describing such behaviors can be exactly solved 
obtaining the complete description of clock states. In addition, the solution 
can be written in an iterative form. 

We therefore provide a clock model including anomalies with an exact 
iterative form that may be useful for simulations as well as for the insertion 
of the clock model in further processing as Kalman hltering, orbit and clock 
estimation algorithms, time scale algorithms, or other types of estimation 
and hltering. Moreover we provide an exact solution that allows the predic¬ 
tion and evaluation of the clock error at a certain time after synchronization 
also in presence of anomalies and as a function of the epoch of the anomaly 
occurrence. This knowledge is of fundamental importance for GNSS evalua¬ 
tions as the clock time deviation directly impact the user solution accuracy 

mn 

In Section 2 we recall the main features of the mathematical model 
based on stochastic differential equations including phase, frequency, and 
drift jumps and we introduce in the model a possible sequence of jumps and 
the increase of instability. In Section 3 we illustrate through examples the 
features captured by the mathematical model. In Section 4 we explicitly ad¬ 
dress the example of a space Rubidium clock on board of a navigation satellite 
and we evaluate the prediction error at a certain time after synchronization 
in presence of a frequency jump. 

To ease the use of this model in simulating clock behavior, we also add 
the availability of a Matlab code that can be freely downloaded from the 
INRIM web site. 
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2 The model of the atomic clock error with 
anomalies 

It is experimentally observed that the signal of an atomic clock can be af¬ 
fected by several anomalies. Our aim is to extend the model proposed in 
introducing anomalies. In the forthcoming we will consider three different 
types of anomalies that involve jumps at deterministic times or changes in 
the variance. 


2.1 Instantaneous jumps 

Let us consider the three-state clock model that we proposed in [7j. The 
three-state clock model X = {(Xi,X 2 ,X 3 )(f),f > 0} is described by the 
three dimensional stochastic differential equation: 

r dXi(f) = (X2(f) + ^i)dt + aidWi{t) 

< dX2{t) = fi2)dt (T2dW2{t) ( 1 ) 

[ dXsit) = usdt + asdWsit) 


with initial condition X(0) = (ci, C 2 , C 3 ), where {Wi{t),t > 0}, z = 1, 2, 3, are 
three indipendent, one-dimensional standard Wiener processes. 

The variable Xi(f) represents the time deviation, the derivative Xi{t) 
represents the frequency deviation of which X 2 {t) is only a component (i.e., 
what is generally called the random walk component). Finally X 3 {t) repre¬ 
sents the frequency drift or aging. 

The metrological interpretation of the model, the solution of ([^ and 
its iterative form can be found in [7]. Here we recall that the probability 
distribution function of the solution X at time t is Normal with mean [7] 
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Let us extend the model of the clock error Q characterized by jumps 
with amplitude a*, i = 1,2,3 in each component: phase, frequency and drift. 
To introduce a simple and treatable mathematical description, we make use 
of general continuous time Markov Chain N(f) = {{Ni, N 2 , N 3 ){t),t > 0}, 
i.e. a continuous-time stochastic process which takes values in some hnite or 
countable set and for which the time spent in each state has an exponential 
distribution (cf. [E]). If we introduce these anomalies into the clock model 
([^ the corresponding mathematical model of the clock becomes 

( dXi{t) = {X2{t) + ni)dt + aidWiit) + aidNiit) 

< dX2^iyt) = (Xslt) + f^2)dt + a2dW2lt) + a2dN3lt) , ( 4 ) 

[ dXs{t) = i^sdt + a^.dW'iit) + a^dNslt) 

where Ni{t), i = 1,2,3 are independent continuous time Markov Chains. 

The easiest particular case considers istantaneous jumps with amplitude 
tti, i = 1, 2, 3 that take place at deterministic times 6i, i = 1,2, 3. It means 
that 

dNi{t) = dH{t — 6i) = 5{t — 6i)dt. (5) 

where H{t) is the Heaviside function and 6{t) is the Dirac delta function 
|18j . Note that the derivatives of the Heaviside function is the Dirac delta 
function, the integral is the so called ramp function and the multiple integrals 
are given by 

r r t” 

/ . . . / H{t) dt. y dt = (6) 

n 

In Figure [^the Heaviside function, its integral and double integral are shown. 

A hrst generalization of this simple model would consider Ni{t), i = 1,2,3 
independent Poisson processes, i.e. anomalies occurring on random epochs, 
the intertimes between the anomalies are exponentially distributed and the 
jumps are deterministic of size 1 HZI. This model will be the subject of a 
future study as it would help in treating the complete clock behaviour for 
example in the analysis of time scale algorithms. 

Here we consider the choice (|^, the system (|^ becomes 

( dXi{t) = {X2{t) + fii)dt + aidWi{t) + aidH{t — 61) 

< dX2{t) = (X^(t) fi2)dt -|- (T2dW2{t) -|- a2dH{t — 62) . ( 7 ) 

[ dX^{t) = ^^dt -|- cT^dW^it) -|- a^dH(t — 6^) 
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Figure 1: Plot of the Heaviside function H{t), its integral and its double 
integral (from below to above). 


Integrating ([^ and recalling 


/ a6{s — fo)ds = aH{t — to) 
Jo 

[ aH{s — to)ds = a{t — to)H{t — to) 
Jo 

[ a{s-to)H{s-to)ds H{t-to), 

Jo I 


( 8 ) 
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the solution can be written in closed form 


Xi(t) — Cl + (c2 + /il)t + (cs + + 0‘lhhl(^) 

+<^2 fo(i — s)dW2(s) + (Js /q 2 '^^ dhF3(s) 

-\-aiH(t — 61 ) + ci 2 (t — 92)H(t — 62 ) 

• ( 9 ) 

X2{t) = C 2 + (cs + /i2)t + + cr2hf2(^) + cr3 /(^(t — s)(ihF3(s) 

+a2-ff(t — ^2) + a3(^ ~ 93)11{t — 93) 

, ^3(^) = C3 + /i3^ + (yzW 3 {t) + a 3 H{t - 6*3) 


The probability distribution function of the solution X at time t is Normal 
with mean M^(t) and covariance matrix S(t), where the mean depends on 
the mean value of the process without anomalies ([^ 


M^(f) = M{t) + 


aiH{t - 9i) + a2it - 92)H(t - 6 ^ 2 ) + 03 ^* H{t - ts) 

(i2H{t — 6 ^ 2 ) + a3{t — 93)H{t — 93 ) 

a3H{t-93) 


( 10 ) 


while its covariance matrix ([^ does not change. 

Let us now consider a hxed time interval [0,T] and an equally spaced 
partition 0 = fo <^i < ••• < t]\f = T and let us denote with r = tk+i — 4, 
k = 0,1,..., iV — 1 the resulting discretization time step. For simplicity 
we hypothesize that the jumps can occur only in the discretization epochs. 
The sampling time can be as small as necessary to ensure this assumption. 
Moreover, observing that 


//( 4+1 -»)- H(h -0)= (11) 

and consequently 

( 4+1 - 9) [Tr(4+i -9)- H{tu - 0 )] = 0 , ( 12 ) 

where is the Kronecker delta [I], we can express the solution at time 

tk+i in terms of the solution at time tk as: 

-^l(^fc+l) = Xi{tk) + (/il + X2{tk))T + (/i2 + X3{tk))^ + Ai3^ + Jk,l 

^ ^ (13 

-^2(^fc+l) = X2{tk) + (/i2 + X3{tk))T + /i3^ + Jk,2 + ®2^(4+i-02) 

X3{tk+i) = X3(tk) + ^3T + Jfc^3 + a3(5(ij,^^_03) 
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where Jfc = (Ja:,i, Jk,2, Jk,z) can be interpreted as an innovation 


ai(iyi(4+i) - W"i(4)) + - W2{h))ds 

+<^3(4+1 —’S)(hF3(s) — W^{tk))ds 

a2(W^2(4+l) - hh2(4)) + - W"3(4))ds 

(X3(W^3(4+i) - Vh3(4)) 


and is Normal distributed with zero mean and covariance matrix Q given by 



ajr + ^2 4 + ^2 4 

^2T + ^3T 

^3 6, 


Q = 


^2^ + ^3T 

^3 2 

(16) 


cr2 4 

*^3 2 

alh 



In the iterative form, we note that the jump acts only once on the corre¬ 
sponding component, then the iterative process carries the effect in the other 
components. 


2.2 Effect of two equal and opposite jumps 

Let us consider the model of the clock error characterized by an anomaly 
on the second component, a temporary frequency jump. The anomaly is 
composed of two istantaneous jumps at times 6q and 6i > 6 q. At time 6q 
we register a positive jump of amplitude a/A and at time 6i we register a 
negative jump of amplitude —a/A, where A = — Oq. The atomic clock 

error can be described by the following stochastic differential equation 


r dX,{t) = {X 2 {t) +fii)dt + a,dW,{t) 

I dX2{t) = {Xsit) + Ii2)dt + a2dW2(t) + ^dHit - 0o) - ^dHit - Oi) (16) 
[ dXslt) = ii^dt + (j^dW'iit). 


Integrating (16) and recalling we get the solution 


Xi(t) — Cl + (c2 + fll)t + (c3 -|- /i2)y + + (^iWiit) 

+(72 I^{t - s)dW2{s) + 173 /o* ^^dW^{s) 

+aH{t — 9i) + ^{t — 6*o)/[6Io,6»i)(^) 

X2{t) = C2 + (c3 -f ^2)t 4/^3^+ (72W2(t) + <73 /o(^ ~ s)dW3{s) 

+ ^Iieo,ei)(t) 

, Asit) = C 3 -A + 3 t + <73143(t) 


(17) 
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where I[a,b)(t) is the window function defined as 


I[a,b){t) = -a)-H{t-b), 


(18) 


and the term +aH{t — 6i) takes care of the phase deviation accumulated 
until time di. 

The probability distribution function of the solution X at time t is Normal 
with mean and covariance matrix S(f). Note that the mean can be 

decomposed in two terms where the hrst corresponds to the mean value of 
the process without anomalies ([^ 


M"(f) = M{t) + 


aH{t - 6 »i) + f (f - 6 »o)/[ 0 o,ei)(t) 
0 


(19) 


while its covariance matrix (j^ does not change. 
The solution in iterative form becomes: 


f ^litk+l) 

S -^2(tfc+l) 

[ ^ 3 ( 4 + 1 ) 


Xl{tk) + (/il + X2{tk))T + (/i2 + X^{tk))\ + P-3W + J] 
X2{tk) + (/^2 + X^{tk))T + P3 V + Jk,2 + X 


Xsitk) + /i3T + 


6 ' 

Sit,. . . -fln'l — S(, 


(tk+l—So) ^{tk+l—Sl) 


.( 20 ) 


where = {Jk,i, Jk, 2 , Jk, 3 ) is (14). Again we note that the effect of the 
jump is present only in the affected component, X 2 in this case, and then 
the iterative process takes care of the effect on the other components. 


Remark 2.1 Note that the algoritms can be easily extended to n determin¬ 
istic jumps. 


2.3 Variance increase 

The atomic clock dynamics can also show changes in the variance. Let us 
model an atomic clock error with a change of the variance values in the 
interval [^o, ^i]. The process for t ^ [^q, ^i] is described as the model without 
anomalies ([^ and its probability distribution is Normal with mean M(f) and 
covariance matrix S(f). At time t G [^ 0 ,^ 1 ] its distribution changes, while 
the mean M(f) is unchanged ([^ its covariance matrix becomes S^^^(f) that 
is the same matrix (|^ but with parameters ai, cr 2 and (T 3 substituted by crj, 
a 2 and cxg. 







This anomaly can be modeled introducing a new innovation matrix in 
the time interval [6o,6i]. The corresponding mathematical iterative model 
becomes 

f -^l(^fc+l) = Xl{tk) + (/^l + X2{tk))T + (/i2 + \ + Jk,l 

j X2{tk+l) = X2{tk) + (/i2 + Xsitk))'^ + h-3 Y '^^.2 ) (21) 

[ -^3(^fc+l) = ^3{tk) + + Jk,3 

where the innovation is 


jfc — {Jk,l, Jk,2j Jk,3) ~ 


iV(0,Q') if4+ie[0o,^i] 
iV(0,Q) if4+i^[0o,^i] 


( 22 ) 


where Q is the covariance matrix (15) while Q' is the same matrix where the 
parameters cxi, a 2 and are substituted by and cig. 


Remark 2.2 In its present form the model accounts for anomalies at deter¬ 
ministic times. In reality anomalies occur at random times. It is possible to 
extend the model and the corresponding algorithms proposed to random times 
9i, i = 0,1,2,3 under the hypothesis that the random variables 6i are inde¬ 
pendent for i = 0,1, 2, 3 and that are independent from the process X without 
anomalies. Under these hypothesis the distribution of the solution can be 
computed conditioning on the time of the anomalies. The iterative solution 
with deteriministic times can still be used for simulations after generating 
random times for the anomalies epochs according to their given distribution. 


3 Simulations and examples 


The simulation of the models introduced in the previous paragraph can be 


easily handled using the iterative solutions (13), (20) and (21). 


In all the cases it is necessary to simulate the innovations that, at any 
instant tk, are independent identical distributed three dimensional Normal 


random variabiles with zero mean and covariance matrix Q given by (15) 


Since the matrix Q is not diagonal, the three components of the vector 
are dependent. Many softwares have routines to generate multidimensional 
Normal random variables. 

To better understand the correlation mechanism, we explicitely describe 
the procedure for generating Jfc in our three dimensional case. It is well 
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known that the vector may be written as = AZ where Z is a standard 
Normal random vector with zero mean and identity covariance matrix I, 
i.e. Z is a vector of three standard Normal independent random variables 
|19j . The matrix A transforms the three components of Z introducing a 
dependency that is captured in the correlation matrix Q, where Q = AA^. 
Here the superscript T denotes the transposed matrix. The values of A can 
be computed using for example Cholesky decomposition method 121 that 
gives 


A 


0 

921 / gii 922 -gii 

V 'Jii 

gSl 932911—931921 

\/<?ll(<?ll<?22-gii) 


0 

0 



932911 — 931921 
\/9ll(9ll922-9ii) . 


(23) 


where i,j = 1, 2,3 are the entry of Q. 

Note that the iterative solutions, useful for simulations, are exact solu¬ 
tions and do not introduce any approximation error. Exact simulations are 
helpful to understand the evolution of the clock behavior or to obtain esti¬ 
mates of functionals of clock states that cannot be computed theoretically. 
The Matlab code can be freely downloaded from the INRIM website. 


3.1 Simulation of anomalous behavior 

Here we report some examples to illustrate the effect of the different types 
of anomalies of the process described in the previous section. 

In Figure]^ the model ([^ with an instantaneous jump on the third com¬ 
ponent is considered and a sample path of the process X(f) is plotted. The 
parameters considered are /ii = /i 2 = /^s = 0, cxi = <72 = <73 = 1 , Oi = 02 = 0, 
03 = 3 and the times of the jump is 63 = 2. The hgure clearly shows the 
jump in the third component, the linear trend of the second component and 
the quadratic trend of the hrst one, as expected by the solution ([^. 

In Figure]^ the model ([^ with instantaneous jumps in all the three com¬ 
ponents is considered and a sample path of the process X(f) is plotted. 
The parameters considered are /ii = /i 2 = = 0, ai = <72 = <73 = 1, 
ai = 02 = 03 = 3 and the times of the jumps are 9i = 6 , 62 = 4, 63 = 2, 
respectively. The hgure clearly shows the jumps in the trajectories and the 
increasing trend of each component sample path after the jumps. 

In Figure]^ the three components of the process X{t) are plotted in case 
of two equal and opposite frequency jumps. The parameters considered are 
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Figure 2: Simulation of the three components of the process X{t) with pa¬ 
rameters /ii = /i 2 = /^3 = 0, (Ti = (72 = ua = 1, ai = 02 = 0, Os = 3 and 
6*3 = 2. The times of the jump is highlighted by an arrow (arbitrary units). 

hi = /72 = h-s = 0, (Ti = (72 = Us = 1, a = 4 and the time interval of the 
jumps is = [4,6]. It is evident the linear trend of the hrst component 

in the interval [4, 6]. 

In Figure]^ the three components of the process X{t) are plotted in case 
of variance increase. The parameters considered are hi = h 2 = ha = 0, 
(Ti = (72 = (73 = 1 , aj = (72 = cTg = 8 , and the time interval of the change 
of variance is [4, 8 ]. The figure shows the effect of an increase of the noise in 
the selected interval. 

The previous hgures show that the model captures the main features of 
the evolution of the atomic clock error. 
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Figure 3: Simulation of the three components of the process X{t) with pa¬ 
rameters /ii = /i 2 = /U 3 = 0, (Ti = ct 2 = CT 3 = 1, ai = 02 = 03 = 3 and 
6*1 = 6 , 62 = 4, 03 = 2. The times of the jumps are highlighted by an arrow 
(arbitrary units). 

4 Application to a space GNSS clock 

In this section we tune the parameters of the model in order to study a case of 
interest concerning a Rubidium Atomic Frequency Standard (RAFS) placed 
on board of a satellite of a Global Navigation Satellite System (GNSS). We 
evaluate the effect of a clock anomaly on the clock prediction error, which 
is a crucial problem in navigation system. As an example we consider the 
standard parameters describing the behavior of a RAFS on board on the 
experimental satellites of the Galileo system IZH. 

The Allan deviation [ 22 ] indicating the stochastic short term noise of 
such RAFS is usually due to white frequency noise and it is assumed to be 
(Ty(r) = 5 ■ The relationship between the Allan deviation and the 
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Figure 4: Simulation of the three components of the process X{t) with pa¬ 
rameters /ii = /r 2 = h-s = 0, cTi = (72 = (73 = 1, a = 4 and the time interval 
of the jumps is [4,6]. The times of the jumps are highlighted by an arrow 
(arbitrary units). 


diffusion coefficient of matrix Q is known and, since 

( 24 ) 

we can infer ai = 5 ■ 10“^^. The other noises are negligible for this example, 
we set (72 = 10“^^, (73 = 10“^^ for numerical reasons only. We assume that 
all the polynomial, linear or quadratic slopes are zero, either because they 
are really negligible or because they have been estimated and appropriately 
compensated. We also assume that the clock is perfectly synchronized at 
epoch zero, therefore all the initial time constants are set to zero. 

With this set of parameter we check which is the time deviation gained 
by such a clock a certain time after synchronization. Note that in the Galileo 
system the clocks are expected to be re-synchronized every 100 min = 6000 
seconds Ell¬ 


is 
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Figure 5: Simulation of the three components of the process X{t) with pa¬ 
rameters = /i 2 = hs = 0 ; (Ti = (72 = cTs = 1 , ( 7 ^ = ( 7 ^ = (Tg = 8 , and the 
time interval of the change of variance is [4, 8 ]. The time interval of variation 
of the variance is highlighted by arrows (arbitrary units). 


Using the clock model without anomalies ([^ the time deviation Xi at 
time t = 6000 s is a Normal random variable with mean value Mi(f) = 
E[Xi{t)] = 0 and standard deviation approximately ^Sn (t) ~ ai\/i = 
5 ■ 10-^^V 60b0 ~ 4 ■ 10-1° s = 0.4 ns. In this case a time deviation interval 
of ±1, 96'y/Sii(f) ~ ±759 ■ lO"!^ s = ±0.8 ns, with 95% conhdence level, 
after 6000 s, represents the prediction error and hence the contribution on 
the positioning error due to unpredictable clock time deviation. Note that 
the polynomial behavior of the three components has not been considered in 
this example. 

To check the effect of a clock anomaly, we consider a frequency jump 
of relative value Frequency jumps are sometimes observed on space 

clock on board of navigation satellites [10]. In our notation a frequency 
jump corresponds to ([^ when we choose 02 = and ai = 03 = 0. We 
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Figure 6: Probability density function and the 95% confidence interval of 
the time deviation Xi at time t = [1000, 3000, 6000,9000] s, with parameters 
= /.'-2 = /Ws = 0, (Ti = 5 ■ 10“^^, cr 2 = CTs = 10“^^ and jump 02 = 10“^^ at 
time 0 = 100 compared with the case without anomalies. 


suppose that the frequency jump occurs at epoch 6* and the effect on the time 
deviation at epoch t depends on the interval {t — 9). Looking to we expect 
the time deviation Xi to be a Normal random variable with mean value 
M^i(f) = a 2 {t — 6) and standard deviation = ai\/i = 5 ■ 

Let us deal with anomaly occurring ai 6 = 100 s or at 0 = 5000 s. In the 
hrst case the anomaly occurs immediately after the synchronization and at 
t = 6000 s its effect is quite impacting. The 95% conhdence level for the time 
deviation Xi gives the interval a 2 (t — 9) ± 1 . 961 ^^ 11 ^ = 5.9 ± 0.8 ns. 

In the second case, the anomaly occurs just at the end of the predictive 
period, and the prediction error at 95% conhdence level in that case is 02 (t — 
9) ± 1.96y^Sii(f) = 1 ± 0.8 ns, therefore the effect is much smaller, very close 
to the case without anomalies discussed above. 
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In Figure is shown the probability density function and the 95% con- 
hdence interval of the time deviation Xi at time t = [1000,3000,6000,9000] 
s, with jump 02 = 10“^^ at time 6 = 100 compared with the case without 
anomalies. 

As we see from (10), the impact of the frequency jump on the time devi¬ 
ation and hence on the time prediction error is linearly depending on t — 6, 
the worst case being at 6^ = 0 giving a conhdence interval a 2 t ± 1.96 

while the best case of no anomaly gives ±1.96^Sii(f) at 95% conhdence 

level. 


5 Conclusion 

A previous mathematical model based on multidimensional stochastic dif¬ 
ferential equation describing the atomic clock error has been extended to 
include the presence of clock anomalies experimentally observed on clocks, 
particularly in GNSS space clocks. The model includes phase, frequency, 
and drift jumps, as well as changes on the random noise power. Examples 
of such trends are presented and the case of a GNSS space clock is discussed 
by evaluating the time deviation gained at a certain time after synchroniza¬ 
tion with or without the effect of a frequency jump. Our results show the 
ability of the proposed model to describe clock anomalies and their effect on 
complex systems. 
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